\documentclass{article}
\usepackage{xltxtra}
\usepackage{ctex}
\usepackage{amssymb}
\usepackage{amsmath}
\usepackage{geometry}
\usepackage{booktabs}
\usepackage{graphicx}
\usepackage{float}
\usepackage{geometry}
\usepackage{booktabs}
\usepackage{tikz}
\usepackage{listings}
\usetikzlibrary{shapes.geometric}
\usetikzlibrary{shapes.arrows}
\usetikzlibrary{positioning}
\usetikzlibrary{graphs}
\usetikzlibrary{quotes}
\usetikzlibrary{arrows}

\makeatother
\geometry{a4paper,top=1.5cm,bottom=1.5cm}
\title{Programming assignments for Section 2}
\author{张皓祥 \\ 3200102536 强基数学2001}
\date{}

\begin{document}
\maketitle

% 构建一个能够实现构造牛顿插值多项式并做运算的项目
\section{项目构建思路}

为了建构一个能够完成构造牛顿插值多项式的库，该项目以仿函数类 \verb|class Function| 和多项式类 \verb|class Polynomial| 
来装载需要接收的信息。具体思路为
\bigskip

\begin{tikzpicture} [node distance = 20pt, bend angle=45, auto]
    \centering
    \node [rounded corners, draw] (pol) {class Polynomial};
    \node [rounded corners, draw, right=35pt of pol]
          (Newton) {class NewtonFormula};
    \node [rounded corners, draw, right=35pt of Newton]
          (fun) {class Function};
    \draw [->] (pol.east) to (Newton.west);
    \draw [->] (fun.west) to (Newton.east);
\end{tikzpicture}

% class Function
\subsection{class Function 仿函数抽象类}
项目中需要涉及到大量函数的运算，除了多项式本身是一类函数，在做插值过程中也会碰到各种函数。为此，本项目中定义了
一个仿函数抽象类。其一方面可以规范所有相关函数的公有特性，另一方面可以为所有仿函数类的传递提供一个统一媒介。
该仿函数类存放于文件 \verb|Function.h| 内，其中仅定义了两个函数：
\begin{verbatim}
    virtual double operator() (const double & _x) const = 0;
    virtual double diff(const double & _x) const;
\end{verbatim}
其中，重载操作符 () 的操作定义为一个纯虚函数，用于规范所有仿函数类都应该对 () 进行重载，并且重载规定输入和输出
样式，实现统一的调用接口。
而函数 \verb|double diff| 则是用于对函数进行求导。该函数作为虚函数可以被覆盖继承，在继承多项式类的时候，该
函数定义即被重新构建。

% class Polynomial
\subsection{class Polynomial 的建构}
由于最终需要实现的程序需要完成生成插值多项式的任务，项目中不可避免地需要对多项式进行一个合理定义。本项目构建
了一个多项式类 \verb|class Polynomial|, 其具有多项式的一般性质，包括构造，求解，线性运算，并且能够打印输出。
该类存放于文件 \verb|Polynomial.h| 中，由于其是一个函数，首先我们对其继承\verb|Function|的虚函数以规范接口。
此外，其主要功能定义如下：

\noindent $\bullet$ 私有成员变量
\lstset{language=C}
\begin{lstlisting}
int deg; // 储存多项式次数
std::vector<double> coff; // 从低次到高次储存多项式系数
\end{lstlisting}
在此，我们把一个多项式的特征定义为次数(deg)和系数(coff)，通过对次数和系数的操作，我们就可以实现多项式
的所有操作。

\noindent $\bullet$ 构造函数
\begin{lstlisting}
Polynomial();   // 初始化构造，把多项式初始化为 0.
Polynomial(const std::vector<double> & vec);
                // 构造一个以 vec 为系数的多项式.
Polynomial(const Polynomial & poly);
                // 复制构造函数，通过深拷贝构造新的多项式类.
\end{lstlisting}

\noindent $\bullet$ 操作符重载
\begin{lstlisting}
double operator() (const double & _x) const;
        // 重载操作符()使其可以对输入值计算对应多项式值。
double operator[] (const int & n) const;
        // 重载操作符[]使其可以取出第n项系数，该操作是只读操作。
void operator= (Polynomial && poly);
        // 左值引用重载操作符 = 把给定 Polynomial 赋值给特定对象。
void operator= (Polynomial & poly);
        // 右值引用重载操作符 = 把给定 Polynomial 赋值给特定对象。
friend Polynomial operator+ (const Polynomial & p1,
                             const Polynomial & p2);
friend Polynomial operator+ (const Polynomial & poly,
                             const double & a);
friend Polynomial operator- (const Polynomial & p1,
                             const Polynomial & p2);
friend Polynomial operator- (const Polynomial & poly,
                             const double & a);
friend Polynomial operator* (const Polynomial & p1,
                             const Polynomial & p2);
friend Polynomial operator* (const double & c,
                             const Polynomial & poly);
// 友元函数重载操作符 +,-,* 使其完成普通多项式的加减乘法
// 以及对常数的运算
\end{lstlisting}
\textbf{附注1：} 在赋值运算 () 过程中，多项式乘法使用的是秦九韶算法，即
$$
a_0+a_1x+a_2x^2+\cdots+a_nx^n = a_0+x(a_1+x(a_2+x(\cdots(a_{n-1}+a_nx)\cdots)))
$$
\textbf{附注2：} 多项式乘法通过数组的乘法操作完成，其数学原理为
$$
\sum\limits_{i=0}^na_ix^i \sum\limits_{j=0}^mb_jx^j = 
\sum\limits_{k=0}^{n+m}c_k x^k, \quad c_k=\sum\limits_{i+j=k}a_i b_j
$$
\noindent $\bullet$ 功能函数
\begin{lstlisting}
double diff(const double & _x) const;
        // 覆盖继承多项式的求导函数运算
Polynomial GetDiff() const;
        // 把求导函数输出
int GetDeg() const;
        // 获取次数
void assign(const std::vector<double> & vec);
        // 按照给定向量 vec 重构该多项式
void PrintCoff() const;
        // 按照系数形式输出多项式的表示
void PrintPol() const;
        // 按照多项式形式输出多项式的表示
\end{lstlisting}
\textbf{附注:} 在求导运算过程中，多项式乘法使用的同样是是秦九韶算法。

\noindent 基于以上的建构，本项目实现了一个占据较小内存完成多项式功能模拟的函数类。

% NewtonFormula
\subsection{class NewtonFormula}
类 \verb|class NewtonFormula| 是用与把给定参数值计算生成牛顿插值公式的类，存放于文件 \verb|NewtonFormula.h| 中。
同时，由于本次作用需要处理特定情况的 Hermite 插值多项式问题，该类也能够计算一类特定情况的 Hermite 插值多项式，
即给定插值节点处的函数值和一阶导数值都为已知条件前提下的 Hermite 插值多项式问题。为了存储数据样本，在这个库内
定义了节点结构体 \verb|struct node|，其中存储了一个节点的自变量值和函数值信息。

\bigskip
\noindent $\bullet$ 私有成员变量
\begin{verbatim}
std::vector<struct node> point; // 使用 node 向量存储所有节点信息
Polynomial poly;                // 根据节点构造的插值多项式
\end{verbatim}
其中\verb|point| 节点用于接收外界节点信息，接收后该类自行计算得到插值多项式。

\bigskip
\noindent $\bullet$ 私有函数
\begin{verbatim}
void ComputePol();               // 生成 Newton 插值多项式
void HermitePol1(const int & N); // 生成 Hermite 插值多项式
\end{verbatim}
\textbf{附注1:} 两种生成插值多项式的基本思路都是计算各节点多项式的差商值，在一个数组上进行计算.
其算法迭代思路为第 k 次迭代，从后往前进行计算
$$
y_{k,i}=\frac{y_{k-1,i}-y_{k-1,i-k}}{x_i-x_{i-k}}, \quad i = k,k+1,\cdots,n
$$
\textbf{附注2:} 两种生成插值多项式的基本思路都是计算各节点多项式的差商值，然后应用秦九韶算法实现
合并同类项，即
\begin{align*}
p(x) & = a_0+a_1(x-x_0)+\cdots+a_n(x-x_0)(x-x_1)\cdots(x-x_{n-1}) \\
& = a_0+(x-x_0)(a_1+(x-x_1)(\cdots(a_{n-1}+a_n(x-x_{n-1}))\cdots))        
\end{align*}

\noindent $\bullet$ 构造函数
\begin{lstlisting}
NewtonFormula(); // 初始化构造函数
NewtonFormula(const std::vector<struct node> & _point);
                 // 根据给定节点值构造插值多项式
NewtonFormula(const Function & func,
              const std::vector<double> & vec);
                 // 根据给定函数和自变量值构造插值多项式
NewtonFormula(const int N, const double * X,
              const double * Y, const double * dY);
                 // 根据给定节点的值构造 Hermite 插值多项式
\end{lstlisting}
\noindent $\bullet$ 公有成员函数
\begin{lstlisting}
double sovle(const double & x) const; // 对x求解插值多项式p(x)
double diff(const double & x) const;  // 对x求解导数值p'(x)
void PrintPoint() const;              // 输出节点信息
void PrintPol() const;                // 打印插值多项式
void PrintPolCoff() const;            // 打印插值多项式系数
Polynomial GetPoly() const;           // 返回插值多项式
Polynomial GetDiff() const;           // 返回插值多项式导函数
\end{lstlisting}
基于此，该库可以实现通过节点和通过函数两种方式建构插值多项式的功能，并且可以计算相应的插值多项式
函数值，以及计算其导数等功能。

\section{项目求解题目应用}

\subsection{B. Test for $f(x)=\frac{1}{1+x^2}$ and the Runge phenomenon.}
本项目中应用写好的牛顿多项式类以完成对函数 $f(x)=\frac{1}{1+x^2}$ 的等距节点插值拟合，分别测试了插值点
个数为 $n = 2,4,6,8$ 这四种情况，插值结果如下：
\begin{figure} [H]
    \centering
    \includegraphics[scale=0.65]{./figure/B_Runge.png} \label{FigB}
\end{figure}
从图中可以看到，随着 $n$ 增大，生成的牛顿插值多项式在 0 附近确实呈现点态收敛的趋势。然而在远离 0 的位置，随着 
$n$ 的增大，函数实际上震荡更加剧烈，如果限制在给定区间上，那么随着 $n$ 逐渐增大，实际上最大差值 
$\Vert p_n - f\Vert _{\infty}=\max\limits_{x\in\left[-5,5\right]}|p_n(x) - f(x)|$ 是递增最终趋向于 
$+\infty$ 的，并且这种震荡仅在靠近端点处出现，从这个趋势就可以发现等距插值拟合函数 $f(x)$ 会出现 Runge 现象。

\subsection{C. Perform Chebyshev interpolation for $f(x)=\frac{1}{1+25x^2}$.}
和取等距插值点不同，本题中选取了 Chebyshev 插值点，即对于给定插值点数 $n$, 选取点 
$x_k=\cos\frac{2k-1}{2n}\pi,k=1,2,\cdots,n$ 作为插值点求相应的插值多项式。对于 $n=5,10,15,20$, Chebyshev 
插值点求解得到的插值多项式图像如下
\begin{figure} [H]
    \centering
    \includegraphics[scale=0.65]{./figure/C_Chebyshev.png} \label{FigC}
\end{figure}
不同于选择等距插值节点的情况，在选择 Chebyshev 节点过程中，随着 $n$ 增大，整体图像是一致收敛于 $f(x)$ 的。
可以得出结论，不同的选择插值节点方式会影响多项式对给定函数收敛的一致性，在选择 Chebyshev 
插值节点时，Runge 现象可以被有效避免。

\subsection{D. Use Hermite interpolation to solve the car problem.}

% Topic (a)
\noindent \textbf{(a). Use a Hermite polynomial to predict the position of the car and its speed for t = 10s.}

把信息 $Time,Distance,Speed$ 存储数值都存入数组内，可调用 \verb|class NewtonFormula| 中构造 Hermite 插值多项式
的构造函数进行运算。

\noindent $\bullet$ 求解得到其 "Distance - Time" 的变化曲线如下：
\begin{figure} [H]
    \centering
    \includegraphics[scale=0.65]{./figure/D_CarPosition.png} \label{FigD1}
\end{figure}
特别地代入 $t=10$ 的条件可以得到对应多项式拟合的值为 $742.503feet$, 根据图像拟合结果也可以得到在 $0-13s$ 内该
车辆近似直线运动。

\noindent $\bullet$ 代入 $t=10$ 求解该多项式对应的导数值，得到速度的近似估计为 $48.3817feet/s$. 这和线性性质
并不相符，究其原因，可以作图如下进行分析：
\begin{figure} [H]
    \centering
    \includegraphics[scale=0.65]{./figure/D_CarSpeed.png} \label{FigD2}
\end{figure}
实际上，导函数是一个18次插值多项式，其线性拟合仅在0-8这一段比较好，在10附近拟合多项式已经出现了震荡，在13以后
导函数就会迅速地偏离且剧烈地70-80这个速度区间，其根本原因还是在于多项式本身具有的震荡性质是不可避免的，因此使用
插值多项式拟合得到的函数其导函数形态拟合程度会比较差。

\bigskip
% Topic (b)
\noindent \textbf{(b). Use the derivation to determine whether the car exceeds 55 mi/h.}

根据图像可以知道，"Speed - Time" 曲线和 "v = 81" 这条直线有多个交点，因此可判断该小车
速度曾到达过 $55mi/h(81feet/s)$. 此外，我们还可以运用牛顿迭代的方法去判断是否存在这样一个时刻。
建立一个库函数，里面定义类 \verb|class NewtonSolver| 作为迭代法求解的求解器，选取5为初值点，可以得到在时刻 
5.91401s 时，对应速度达到了 $81feet/s$.

\subsection{E. Newton interpolation for tree problem.}

% Topic (a)
\noindent \textbf{(a). Approximate the average weigh curve for each sample.}

分别把两棵树的 “Weight - Day” 关系代入进行牛顿插值运算，绘图可得其关系为如下：
\begin{figure} [H]
    \centering
    \includegraphics[scale=0.65]{./figure/E_Tree29days.png} \label{FigE1}
\end{figure}
如图，即为使用插值多项式得到的曲线模拟图。注意到模拟函数在两端出现了极大的震荡偏差，因此其中合理的模拟区间可以认为
仅有 [0,29] 这一段。

\bigskip
% Topic (b)
\noindent \textbf{(b). Predict whether the two samples will die after another 15 days.}

首先使用模拟的插值多项式来观察，如图
\begin{figure} [H]
    \centering
    \includegraphics[scale=0.65]{./figure/E_Tree44days.png} \label{FigE2}
\end{figure}
图像在30天以后对曲线的模拟已经出现了严重偏移，不能使用多项式图像作为判断依据了。为了更加准确的进行预测，我们应该
仅仅选取0-29天的函数曲线作为判断依据。而在这一段时间，两棵树的近似曲线都在10-29这一段内呈现下降趋势。沿着这个趋势
进一步作图，可以预测其会在 $x<44$ 时出现零点，因此可以合理预测在另一个15天以后两棵树都死亡了。

\section{总结}

该项目实现了对给定数据生成牛顿插值多项式的功能，并实现了计算特定条件下的 Hermite 插值多项式的功能。通过该项目，
我们探究了 Runge 现象的出现问题，并且通过合理选取 Chebyshev 插值节点有效规避了 Runge 现象的出现。接下来，该
项目有望进一步改进优化，比如完成对一般的 Hermite 插值多项式的计算，以及探讨更多相关的插值问题。

\end{document}